#include "includesDPC.h"
#include "pocsDPC.h"
// Code performing Hilbert projection.
using namespace std;
void P1hilbert(ScanData & s)
{
	int i,j,n;
	double pterm=0,x=0,a = 0.60 ,b = 1.-a,temp2=0,temp=0;
	vector<double> image (s.Ny, 0);
	for (j=0;j<s.Ny;j++){
		// Compute preterm
		for (i=0; i<s.Nx;i++){
			pterm = pterm + s.dx*s.recon[i][j];
		}
		//cout<< pterm<<endl;
		// Compute H transform
		for(i=0; i<s.Nx;++i){ // x loop
			image[i]=0;
			for(n=0; n<s.Nx;++n){ // x prime loop
				if( ((i-n)%2 != 0) ) image[i] = image[i] + s.recon[n][j] * (2./(n-i))/M_PI;
			}
			s.HT[i][j] = image[i];
		}
			// perform smooting

		temp = image[0];
		for(i=1; i<s.Nx-1;++i){ // x loop
			temp2=image[i];
			image[i] = 0.5*b*temp + a*image[i] + 0.5*b*image[i+1];
			temp=temp2;
		}

		
		// Compute projection

		for(i=0; i<s.Nx;++i){
			if (s.DBPmask[i][j]) image[i] = s.DBPdata[i][j]*(-0.5/M_PI) + 0.02 ;// + 0.04354751;
		}
		
		// Compute inverse H transform
		for(i=0; i<s.Nx;++i){ // preweight
			x = s.xmin + i*s.dx;
			image[i] = image[i] * sqrt(s.X2 - x*x);
		}
		
		for(i=0; i<s.Nx;++i){ // convolution
			s.recon[i][j]=0;
			for(n=0; n<s.Nx;++n){ // x prime loop
				if( ((i-n)%2 != 0)) s.recon[i][j] = s.recon[i][j] + image[n] * (-2./(n-i));
			}
			
		}

		// perform smooting
		temp = s.recon[0][j];
		for(i=1; i<s.Nx-1;++i){ // x loop
			temp2=s.recon[i][j];
			s.recon[i][j] = 0.5*b*temp + a*s.recon[i][j] + 0.5*b*s.recon[i+1][j];
			temp=temp2;
		}

		for(i=0; i<s.Nx;++i){ // postweight
			x = s.xmin + i*s.dx;
			s.recon[i][j] = (pterm + s.recon[i][j])/sqrt(s.X2 - x*x)/M_PI;
			//s.recon[i][j] = (s.recon[i][j])/sqrt(s.X2 - x*x)/M_PI;
			image[i]=0;
		}
		
		pterm =0;
	}
}
